function [theta_new, lag_new] = search_for_theta(model_sol, lower_spread_target)
model_sol.par.model_name='behavioral';  % Must be behavioral to make sure that it matters

theta_upper = 5;  theta_lower=0; 
theta_new = model_sol.par.theta_diagnostic;
theta_last = theta_new;  diff=1;

iter=1;  iter_count_max = 10;
while( abs(diff) > 0.01 && iter <iter_count_max )
    m = simulation_results_output(model_sol,8000); 
    coef= m.reg_spread_too_low_with_control.Coefficients.Estimate(2,1)/std(m.credit_spread);
    if(  coef<=lower_spread_target )
        theta_upper = model_sol.par.theta_diagnostic;
    else
        theta_lower = model_sol.par.theta_diagnostic;
    end
    theta_new = (theta_upper+theta_lower)/2;
    theta_last = model_sol.par.theta_diagnostic;
    diff = abs(coef-lower_spread_target)/abs(lower_spread_target);
    disp( [ 'Credit spread too low by ', num2str(coef), ' with theta=',  num2str(model_sol.par.theta_diagnostic)] );
    iter = iter+1;
    model_sol.par.theta_diagnostic = theta_new;
end

lag_new=1;
if( abs(coef-lower_spread_target)/abs(lower_spread_target) > 0.1  )
    lag_upper = 5;  lag_lower=0.5; 
    lag_new = model_sol.par.lag_years;
    lag_last = lag_new;  diff=1;
    iter=1;  iter_count_max = 10;
    while( abs(diff) > 0.01 && iter <iter_count_max )
        m = simulation_results_output(model_sol, 8000); 
        coef= m.reg_spread_too_low_with_control.Coefficients.Estimate(2,1)/std(m.credit_spread);
        if(  coef<=lower_spread_target )
            lag_upper = model_sol.par.lag_years;
        else
            lag_lower = model_sol.par.lag_years;
        end
        lag_new = (lag_upper+lag_lower)/2;
        lag_last = model_sol.par.lag_years;
        diff = abs(coef-lower_spread_target)/abs(lower_spread_target);
        disp( [ 'Credit spread too low by ', num2str(coef), ' with lag=',  num2str(model_sol.par.theta_diagnostic)] );
        iter = iter+1;
        model_sol.par.lag_years = lag_new;
    end
end

end